/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#if !defined _Cauchy_General_Stochastic_Tree_
#define _Cauchy_General_Stochastic_Tree_

#include <src/math/Stochastic/Stochastic_Tree.h>
#include <src/algorithm/Evolution_Algorithm/all.h>

/*! \brief stochastic processes transport map solved through a Crank Nicolson algorithm.

Consider a stochastic process driven by the following SDE
\f[
    dX = r(t,X) dt + \sigma(t,X) dW
\f]
with \f$X(0)=S_0\f$ and \f$X\in \Omega\f$. Denote \f$\rho(t,x)\f$ its density, \f$S^{-1}(t,x) = \int_{s\in\Omega;s<x} \rho(t,s)ds\f$ its cumulative
and denotes \f$S(t,y)\f$, from \f$[0,1]\f$ into \f$\Omega\f$ its inverse map. \f$S\f$ follows the PDE
\f[
    \partial_t S = r(t,S) - \frac{1}{2} \partial_y \left(\frac{\sigma^2(S)}{\partial_y S} \right)
\f]
Since the density \f$\rho(t,S) = \rho(t,S(y))\f$ follows
\f[
    \partial_t \rho(S) = \rho^2(S) \left( \frac{1}{2} \partial_{y} \left( \sigma^2(S) \partial_y \rho(S)\right)  - \partial_y \left(r(S) - \frac{1}{2} (\sigma^2)'(S)\right) \right)
\f]
This class proposes a Crank Nicolson type algorithm to solve the previous equation.
Let us introduce some notations. Consider the discrete space grid
\f$\{y_i = \frac{i+1/2}{N}\}_{i=1,...,N-1}\f$ and a discrete time grid \f$\{t^i\}_{i=1,..,D}\f$ (this set will be made precise below).
Denote \f$\delta_{\bar{y}}\f$ (resp. \f$\delta_{y}, \delta_{\hat{y}}\f$) the backward
(resp.forward, centered) finite difference operator \f$\delta_{\bar{y}}f_i =
\frac{f_{i}-f_{i-1}}{h}=N\left(f_{i}-f_{i-1}\right)\f$ (resp. \f$\delta_{y}f_i =N\left(f_{i+1}-f_{i}\right)\f$, \f$\delta_{\hat{y}}f_i =2N\left(f_{i+1}-f_{i-1}\right)\f$). The discrete Laplacian is denoted \f$\delta_{y \bar{y}}^2 = \delta_{y}\delta_{\bar{y}}\f$.
The forward time derivative is denoted \f$\delta_t f^n = \frac{f^{n+1}-f^n}{t^{n+1}-t^n} := \frac{f^{n+1}-f^n}{\tau^n}\f$.
We denote the discrete map \f$S_i^n \equiv S(t^n,y_i)\f$, and \f$S_i^{n+1/2} = \frac{S_i^n+S_i^{n+1}}{2}\f$. For any function \f$f(S)\f$,
we denote \f$f_i^n := f(S_i^n)\f$ and \f$f_i^{n+1/2} = f \left(S_i^{n+1/2} \right)\f$.

Let us consider the following scheme:
\f[ \delta_t S_i^n = \tilde{r}_i^{n+1/2} - \frac{ \left( \sigma_{i}^{n+1/2} \right)^2}{2} \delta_{\hat{y}} \left( \rho_i^{n+1/2} \right) \f]
\f[ \rho_i^{n+1/2} = \frac{1}{\delta_{\hat{y}} S_i^{n+1/2}} \f]

This scheme in unconditionally stable since we can prove that it verifies the following discrete entropy evolution :
    \f[
        \delta_t \|S^n\|_{\ell^2}^2 = \|\sigma^{n+1/2}\|_{\ell^2}^2 +2 \tilde{r}^{n+1/2} \cdot  S^{n+1/2}
    \f]
where we denoted \f$\|S^{n}\|_{\ell^2}^2 = \frac{1}{N}\sum_i (S_i^n)^2\f$, and \f$\frac{1}{N} f\cdot g = \sum_i g_i f_i\f$.
This last inequality provides the bound \f$\|S^n\|^2 \le C\exp(t^n)\|S^0\|^2\f$ under the sub-quadratic hypothesis (\ref{drift}).
This non linear scheme may be solved by a generic iterative process : let us define \f$S^{n+1,0} = S^n\f$, then compute recursively for \f$k\ge 0\f$
    \f[ S_i^{n+1/2,k-1} = \frac{S_i^n+S_i^{n+1,k-1}}{2} \f]
    \f[ \rho_i^{n+1/2,k-1} = \frac{1}{\delta_{\hat{y}} S_i^{n+1/2,k-1} }  \f]
    \f[  S_i^{n+1,k} = S_i^{n} + \tau^n \left( \tilde{r}_i^{n+1/2,k-1} - \frac{\left(\sigma_i^{n+1/2,k-1}\right)^2}{2}\delta_{\hat{y}}\rho_i^{n+1/2,k-1} \right) \f]
We end the procedure until a criteria is reached, namely
\f[
    \frac{\| S^{n,k+1}- S^{n,k}\|_{\ell_2}}{\|S^{n}\|_{\ell_2}} \le 0.01
\f]
Note that, even if this scheme is unconditionally stable, the iterative algorithm used is convergent under a C.F.L. condition :
\f[
    \tau^n N^2 \le C \| \rho^n (\sigma^n)^2\|_{\ell^{\infty}}
\f]
*/
	template <class stochastic_process>
	class Crank_Nicolson_Sto : public Crank_Nicolson <
		Delta_t<typename stochastic_process::underlying_vector_type, 
		stochastic_process::underlying_vector_type::Dim>>
	{
	public:
		typedef typename stochastic_process										stochastic_process_type;
		typedef typename stochastic_process_type::underlying_type				underlying_type;
		typedef typename stochastic_process_type::underlying_vector_type		underlying_vector_type;
		typedef typename stochastic_process_type::density_vector_type			density_vector_type;

		typedef typename OS_Multi_Array<Neumann_Vector,
										stochastic_process_type::Dim+1,
										typename stochastic_process_type::base_value_type>::type
														rate_vector_type;

		typedef Crank_Nicolson Base;

//!compute the time steps as 
//!	\f[
//!    \tau^n =\frac{C^n}{N^2} \| \rho^n (\sigma^n)^2\|_{\ell^{\infty}}
//!\f]
//! where C^n is given by a time adaptive technique.
		virtual OS_double CFL() {
			densities_  = 1./u_c_(S_k12);
			norm2_ = jmmath::norm22(iterate);
			diffusion_ = process_->sigma_helper<underlying_vector_type,underlying_vector_type>(t_12,iterate);
			double CFL = jmmath::norm_infty(diffusion_^densities_^densities_);
			double size = S_k12.size();
			CFL = 1./CFL/size/size;
			delta_t_ =  std::min(delta_t_, cfl_adjust_*CFL);
/*			double drift = std::exp(process_->mu(current_time_, process_->x0())*delta_t_);
			OS_Array_pow(iterate, std::sqrt(std::sqrt( 1. + delta_t_/current_time_) ) );
			OS_Array_multiply(iterate,drift);*/
			return CFL;
		};
//! return true if
//!
//!		We end the procedure until a criteria is reached, namely
//! \f[
//!    \frac{\| S^{n,k+1}- S^{n,k}\|_{\ell_2}}{\|S^{n}\|_{\ell_2}} \le 0.01
//! \f]
		virtual bool criteria() {
			OS_double breakpoint1 = jmmath::norm22(iterate - S_k1);
//			OS_double breakpoint1 = jmmath::norm22(Op_S_k1 - Op_iterate);
			return breakpoint1/(1.+norm2_) < 1e-5;
		};
//!    computes \f[\rho_i^{n+1/2,k-1} = \frac{1}{\delta_{\hat{y}} S_i^{n+1/2,k-1} } \f] and \left( \tilde{r}_i^{n+1/2,k-1} - \frac{\left(\sigma_i^{n+1/2,k-1}\right)^2}{2}\delta_{\hat{y}}\rho_i^{n+1/2,k-1} \right)]
		virtual void gstar() {
			diffusion_ = process_->sigma_helper<underlying_vector_type,underlying_vector_type>(t_12,S_k12);
			rate_ = process_->mu_helper<underlying_vector_type,rate_vector_type>(t_12,S_k12);
			densities_12  = 1./u_c_(S_k12);
			rhosigma_12 = diffusion_^densities_12;
			Op_iterate = rate_-d_c_(rhosigma_12);
		};
//!		compute \f[S_i^{n+1/2,k-1} = \frac{S_i^n+S_i^{n+1,k-1}}{2}\f]
		virtual void vstar() {
			S_k12   = (*get_state_variable()+iterate)/2.;
		};
		void set_process(smart_ptr<stochastic_process_type> & process) {
			process_ = process;
		};
		virtual void update(){
			Base::update();
		};
	virtual void adjust_cfl(const time_type & from, argument_type S_k){
			Base::adjust_cfl(from,S_k);
/*			double drift = std::exp(process_->mu(current_time_, process_->x0())*delta_t_);
			OS_Array_pow(iterate, std::sqrt(std::sqrt( 1. + delta_t_/current_time_) ) );
			OS_Array_multiply(iterate,drift);*/
	}
		smart_ptr<stochastic_process_type>	process_;
		underlying_vector_type				diffusion_;
		rate_vector_type					rate_;
		density_vector_type					densities_;
		density_vector_type					densities_0;
		density_vector_type					rhosigma_12;
		density_vector_type					densities_12;
		OS_double							norm2_;
		D_c_<density_vector_type>			d_c_;
		D_c_<underlying_vector_type>		u_c_;
	};

/*! \brief Cauchy Interface to handle the Cauchy problem linked to Crank_Nicolson_Sto.
*/	
	template <class tree_type>
	class Cauchy_Stochastique : public CauchyProblem_traits< 
		Crank_Nicolson_Sto<tree_type> , 
		Graphical_Interface>
	{
		typedef typename tree_type::stochastic_process_type stochastic_process_type;
		public :
			Cauchy_Stochastique(smart_ptr<stochastic_process_type> &process) : process_(process)
			{
				get_algorithm().set_process(process_);
			};
			virtual void set_dim(OS_size dim) {
				OS_DYNAMIC_CHECK(dim == 1,"Cauchy_StochastiqueOnly Dim 1");
			};
			OS_size get_dim() {return 1;};
			void initial_conditions(OS_size nb_part)
			{
			};
			virtual smart_ptr<interface_type> Create_Instance() {
				return smart_ptr<interface_type>(new Cauchy_Stochastique(process_));
			};
		smart_ptr<stochastic_process_type> process_;
	};



#endif